#include <algorithm>
#include "manager.h"
#include "aggAux.h"

// Performs the Super Edge creation and coefficients addition
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// IN: theMesh, setOfDomains, iMesh
// OUT: aggMesh

void treatBoundary (pMesh theMesh, pEdge edge, pEdge superEdge, int seedID_1, int seedID_2, pVertex node_1, pVertex node_2,
                    set<int>::iterator iter);

void aggSE(pMesh theMesh, set<int> &setOfDomains, vector<pMesh> aggMesh, int iMesh)
{

cout<<"   Creating Super Edges... ";

// Definig the domain flag
set<int>::iterator iter = setOfDomains.begin();

int iSE(0); // Superedge counter

pMesh fineMesh = (iMesh == 0)?theMesh:aggMesh.at(iMesh - 1);

// Variables used to store the C's coefficients
double val, val2;
char *Cij_string[2];
for (int i=0; i<2; i++)
{
	Cij_string[i] = new char[64];
}

// Loop on the fine mesh edges
EIter eit = M_edgeIter(fineMesh);
while (pEntity edge = EIter_next(eit))
{
	// Get edge´s node
	pVertex node_1 = E_vertex(edge,0);
	pVertex node_2 = E_vertex(edge,1);

	// Get node´s seeds IDs
	int seedID_1(-1), seedID_2(-1);
	EN_getDataInt(node_1, MD_lookupMeshDataId("SeedID"), &seedID_1);
	EN_getDataInt(node_2, MD_lookupMeshDataId("SeedID"), &seedID_2);

	if (seedID_1 == seedID_2)
	{   // Do nothing! Same seed for both nodes, so that this edge will not exist
		// in the coarse mesh.
	}else
	{   // Different seeds! So treat this edge will exist in the coarse mesh
		// Ordering it so that (node_1(ID)<node_2(ID))
		int minID(0), maxID(0);
        float signal(0);

		if (seedID_1 < seedID_2)
		{
		    minID = seedID_1;
		    maxID = seedID_2;
		    if (EN_id(node_1) < EN_id(node_2))
		    {
                signal = 1.;
		    }else
		    {
		        signal = -1.;
		    }
		}else
		{
		    minID = seedID_2;
		    maxID = seedID_1;
		    if (EN_id(node_2) < EN_id(node_1))
		    {
                signal = 1.;
		    }else
		    {
		        signal = -1.;
		    }
		}

		// Get node´s seeds vertexs
		// (DEV) - VERIFICAR A DIFERENCA ENTRE mVERTEX
		// e pVERTEX!!!
		pVertex minSeed = (pVertex)aggMesh.at(iMesh)->getVertex(minID);
		pVertex maxSeed = (pVertex)aggMesh.at(iMesh)->getVertex(maxID);

        if (aggMesh.at(iMesh)->getEdge((mVertex*)minSeed, (mVertex*)maxSeed) == 0)
		{   // Edge does not exist yet

			pEdge newSuperEdge = M_createE(aggMesh.at(iMesh), minSeed, maxSeed, 0);
			// Associate Cij coefficient to edge
			// (DEV) - VERIFICAR COMO OBTER gEntity!!
			sprintf(Cij_string[0],"Cx-%d",*iter);
            sprintf(Cij_string[1],"Cy-%d",*iter);
			for (int j=0; j<2; j++)
			{
				val = .0;
				EN_getDataDbl(edge,MD_lookupMeshDataId(Cij_string[j]),&val);
				val *= signal;
				EN_attachDataDbl(newSuperEdge,MD_lookupMeshDataId(Cij_string[j]),val);
			}

			// Determine if the edge is on geometry or domain boundary,
			// in order to add the D coefficients
			// Looking for boundary flag
			int isBoundary(0);
			EN_getDataInt (edge, MD_lookupMeshDataId("isBoundary"), &isBoundary);

			if (isBoundary)
			{
			    treatBoundary(fineMesh, edge, newSuperEdge, seedID_1, seedID_2, node_1, node_2, iter);
			}
			iSE++;
		}else
		{   // Edge already exists

			pEdge oldSuperEdge = (pEdge)aggMesh.at(iMesh)->getEdge((mVertex*)minSeed, (mVertex*)maxSeed);

			// Associate Cij coefficient to edge
			// (DEV) - VERIFICAR COMO OBTER gEntity!!
			sprintf(Cij_string[0],"Cx-%d",*iter);
            sprintf(Cij_string[1],"Cy-%d",*iter);

			for (int j=0; j<2; j++)
			{
				val = .0;
				val2 = .0;
				EN_getDataDbl(edge,MD_lookupMeshDataId(Cij_string[j]),&val);
				val*= signal;
				EN_getDataDbl(oldSuperEdge,MD_lookupMeshDataId(Cij_string[j]),&val2);
				val2 += val;
				EN_attachDataDbl(oldSuperEdge,MD_lookupMeshDataId(Cij_string[j]),val2);
			}

			// Looking for boundary flag
			int isBoundary(0);
			EN_getDataInt (edge, MD_lookupMeshDataId("isBoundary"), &isBoundary);

			if (isBoundary)
			{
                treatBoundary(fineMesh, edge, oldSuperEdge, seedID_1, seedID_2, node_1, node_2, iter);
			}

		}
	}
}
EIter_delete(eit);
cout<< iSE << " edges in the "<< iMesh <<"° coarse mesh \n";
}

// Routine to compute the D coefficients of the boundary edges
void treatBoundary (pMesh fineMesh, pEdge edge, pEdge superEdge, int seedID_1, int seedID_2, pVertex node_1, pVertex node_2,
                    set<int>::iterator iter)
{

    // Variables used to store the coefficients
    double val, val2, val3, valTotal, deltaDij;

    char *Dij_string[2], *deltaDij_string[2];

    for (int i=0; i<2; i++)
{
	Dij_string[i] = new char[64];
	deltaDij_string[i] = new char[64];
}

    // Attaching boundary flag
	EN_attachDataInt(superEdge, MD_lookupMeshDataId("isBoundary"), 1);

	// Get node´s seeds vertex
    pVertex seed_1 = (pVertex)fineMesh->getVertex(seedID_1);
	pVertex seed_2 = (pVertex)fineMesh->getVertex(seedID_2);

	pEdge boundEdge1 = (pEdge)fineMesh->getEdge((mVertex*)node_1, (mVertex*)seed_1);
	pEdge boundEdge2 = (pEdge)fineMesh->getEdge((mVertex*)node_2, (mVertex*)seed_2);

    // Associate D coefficients to edge
	// (DEV) - VERIFICAR COMO OBTER gEntity!!
	sprintf(Dij_string[0],"Dx-%d",*iter);
    sprintf(Dij_string[1],"Dy-%d",*iter);
    sprintf(deltaDij_string[0],"deltaDx-%d",*iter);
    sprintf(deltaDij_string[1],"deltaDy-%d",*iter);
    for (int j=0; j<2; j++)
	{
		val = .0;
		val2 = .0;
		val3 = .0;
		valTotal = .0;

        if (boundEdge1 != 0)
		{
            deltaDij = .0;
			EN_getDataDbl(boundEdge1, MD_lookupMeshDataId(Dij_string[j]),&val);
			EN_getDataDbl(boundEdge1, MD_lookupMeshDataId(deltaDij_string[j]),&deltaDij);
        }
		if (boundEdge2 != 0)
		{
            deltaDij = .0;
			EN_getDataDbl(boundEdge2, MD_lookupMeshDataId(Dij_string[j]),&val2);
			    EN_getDataDbl(boundEdge2, MD_lookupMeshDataId(deltaDij_string[j]),&deltaDij);
					}

                    deltaDij = .0;
					EN_getDataDbl(edge, MD_lookupMeshDataId(Dij_string[j]),&val3);
					EN_getDataDbl(edge, MD_lookupMeshDataId(deltaDij_string[j]),&deltaDij);
					double val3MD = val3; // Coefficient according to Median Dual
					if (seedID_1 < seedID_2)
					{
                        if (EN_id(node_1) < EN_id(node_2))
                        {
                            val3 += deltaDij;
                        }else
                        {
                            val3 -= deltaDij;
                        }
                    }else{
                        if (EN_id(node_2) < EN_id(node_1))
                        {
                            val3 += deltaDij;
                        }else
                        {
                            val3 -= deltaDij;
                        }
                    }

                    // D coefficient according to median dual (Dij = Dji)
                    valTotal = val + val2 + val3MD;

                    // Calculate delta D to correct the error that arises when it is assumed that the
                    // coarse mesh is also median dual, what is not necessarily true (Dij != Dji)
                    double newDeltaDij = .0;
                    if (seedID_1 < seedID_2)
                    {
                        newDeltaDij = 2*val + val3 - valTotal;
                    }
                    else // minID == seedID_2
                    {
                        newDeltaDij = 2*val2 + val3 - valTotal;
                    }

                    EN_attachDataDbl(superEdge,MD_lookupMeshDataId(deltaDij_string[j]),newDeltaDij);
                    EN_attachDataDbl(superEdge,MD_lookupMeshDataId(Dij_string[j]),valTotal);
				}
}
